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Foreword 

This  report  was  prepared  by  Professor  Richard  H.  Rapp,  and  Jaime  Y. 
Cruz,  Post-Doctoral  Researcher,  Department  of  Geodetic  Science  and  Surveying, 
The  Ohio  State  University,  under  Air  Force  Contract  No.  F19628-86-K-0016,  OSU 
Project  No.  718188.  The  contract  covering  this  research  is  administered  by 
the  Air  Force  Geophysics  Laboratory,  Hanscom  Air  Force  Base,  Massachusetts, 
with  Dr.  Christopher  Jekeli,  Scientific  Program  Officer. 

The  computer  resources  provided  at  the  Boeing  Computer  Center  were 
provided  by  a  grant  from  the  National  Science  Foundation,  EAR-8420862. 

Computer  resources  at  The  Ohio  State  University  were  provided  by 
contract  funds  and  by  the  Instruction  and  Research  Computer  Center. 

The  orbit  calculations  were  carried  out  at  NASA's  Goddard  Space  Flight 
Center  with  the  cooperation  and  help  of  Mr.  James  Marsh  and  Mr.  Steve 
Klc  ?ko. 

Mr.  Vasilios  Despotakis  was  most  helpful  in  carrying  out  the  estimation  of 
the  l*xl*  terrestrial  mean  anomaly  field. 


1.  In  troductio  n 

The  development  of  high  degree  spherical  harmonic  fields  has  been  carried 
out  for  a  number  of  years.  The  estimation  of  such  fields  has  become  more 
viable  due  to  the  availability  of  satellite  altimeter  data  in  the  ocean  areas,  and 
the  increased  availability  of  terrestrial  gravity  measurements.  One  of  the  first 
solutions  to  degree  180  was  that  of  Rapp  (1978).  In  1981  Rapp  carried  out  a 
second  expansion  to  degree  180  using  a  more  complete  set  of  a  priori  potential 
coefficients  than  used  in  the  1978  solution.  Other  solutions  have  been 
developed  by  Lerch  et  al  (1981,  GEM10C)  and  Wenzel  (1985,  GPM2). 


In  the  past  few  years  our  studies  have  been  related  to  improved  modeling 
techniques  taking  into  account  approximations  made  earlier  (Rapp,,  1984;  Cruz, 
1985)  and  in  the  improvement  of  our  l'xl*  mean  anomalies  in  land  and  ocean 
areas  (Rapp,  19P6a).  The  progress  in  these  areas  has  led  us  to  the 
development  of  new  high  degree  expansions  that  take  advantage  of  the  new 
theory  and  new  data.  This  report  describes  the  data,  methods,  and  analysis 
used  in  the  development  of  the  new  solutions. 


2±  Data  to  be  Used 

We  first  describe  each  of  the  possible  sources  of  data  for  use  in  the  new 

sol  u  t  ion  s. 


2.1  l'xl*  Terrestrial  Gravity  .Data 

W.«  hav--  r<',-rnt!y  completed  an  update  of  our  l'xl'  mean  free  air  anomaly 
data  set.  (Despotnkis,  1986).  This  data  set  has  evolved  from  prior  compilations 
and  from  n  merger  with  data  made  available  from  DMA  Aerospace  Center.  The 
new  data  -et  contains  48,955  l'xl'  anomalies  (Figure  1)  and  5689  anomalies 
estimated  by  geophysical  correlation  techniques.  The  location  of  most  of  the 
geophysical  anomalies  is  shown  in  Figure  5.  This  new  data  set  has  corrected 
a  number  of  e-rors  found  in  the  set  used  for  the  OSU81  solution  and  has 
incorp'  ra'ed  »  number  of  new  date  sources.  The  new  data  has  led  to  better 
anomaly  estimates  for  much  of  Europe,  all  of  Australia,  parts  of  Greenland, 
South  Africa.  Tnd  a.  Japan,  Brazil,  and  other  areas. 


2.2  Altimeter  Derived  l'xl*  Gravity  Anomalies 

A  combined  Geos-3/Seasat  altimeter  data  set  has  been  used  to  derive 
gravity  anomalies  on  n  0'.125  grid  in  the  ocean  areas.  This  gridded  data  has 
been  used  to  -■  Unii  30  x90'  and  l'xl'  mean  free  air  anomalies  (Rapp,  1986a). 
Approximately  37,000  values  were  estimated  although  some  values  were  on  land 
and  some  values  were  not  reliable  because  of  the  sparsity  of  altimeter  data  in 
the  area.  The  location  of  all  anomalies  is  shown  in  Figure  2.  In  later  use  we 
will  d e ! e * <>  fron  tj,.s  data  set  anomalies  on  land  and  values  with  large  standard 
deviation*--. 


2.  Use  a  terrestrial  anomaly  in  ocean  areas  only  if  its  standard  deviation  is 
<  5  mgals  or  if  no  altimeter  derived  anomaly  exists. 

3.  Do  not  use  altimeter  derived  anomalies  south  of  latitude  -65*. 

4.  For  accepted  terrestrial  data  apply  the  following  corrections: 

A.  Due  to  Mass  of  Atmosphere  (see  equation  (3.43)) 

B.  To  Reduce  Surface  Anomaly  to  the  Ellipsoid  (see  equation  (3.30)) 

C.  For  Gravity  Formula  Change 

5.  Reduce  altimeter  derived  anomalies  to  the  ellipsoid  from  the  geoid  using  a 
gradient  technique  (see  (3.30)). 

6.  Do  not  use  altimeter  derived  anomalies  in  the  Mediterranean  area  defined 
as  (30*  <  ♦  *  50*;  0*  *  X  *  40*). 

7.  Create  two  data  sets  where  one  excludes  geophysical  anomalies  and  the 
other  includes  such  anomalies.  Preliminary  tests  indicated  that  33  geophysical 
anomalies  were  very  reasonable  values  and  they  were  therefore  included  in 
the  first  data  file.  The  location  of  these  anomalies  and  the  rational  for  their 
selection  will  be  discussed  later.  We  will  be  making  two  main  combination 
solutions;  one  will  be  designated  OSU86C  and  it  will  exclude  geophysical 
anomalies  (except  for  the  select  33  values);  the  other  will  include  the 
geophysical  anomalies  and  will  be  designated  OSU86D. 

The  location  of  the  50,562  l*xl*  anomalies  used  for  the  OSU86C  solution  is 
shown  in  Figure  3  while  Figure  4  shows  the  location  of  the  56,109  values  used 
in  the  OSU86D  solution.  Figure  5  shows  the  location  of  5,547  l*xl*  anomalies 
used  in  OSU86D  but  not  in  OSU86C.  Figure  6  shows  the  location  of  the  32,274 
anomalies  derived  from  the  altimeter  data  that  are  used  in  both  the  C  and  D 
solutions. 

In  order  to  form  a  global  field,  values  have  to  be  estimated  for  blocks  in 
which  r.o  anomaly  estimate  is  available.  Setting  such  values  to  zero  is 
inconsistent  with  knowledge  of  the  gravity  field  that  we  do  have  from  the  a 
priori  potential  coefficient  estimates.  We  therefore  calculated  a  rigorous 
anomaly  on  the  ellipsoid  using  the  potential  coefficients  described  in  Section 
2.3. 


The  end  product  of  this  merger  is  two  sets  of  64,800  l*xl*  anomalies  that 
refer  to  an  ellipsoid  whose  parameters  were  defined  earlier  and  such  that 
there  is  no  atmosphere  outside  the  ellipsoid.  Each  anomaly  has  been  assigned 
a  standard  deviation  which  comes  from  the  original  estimation  process.  In  the 
case  of  the  anomaly  "fill  ins”  we  assigned  a  standard  deviation  of  *25  mgal. 
It  is  our  intention  to  assume  that  all  the  anomaly  estimates  are  independent. 
This  is  not  the  case  and  problems  caused  by  our  assumption  will  be  discussed 
later. 

We  give  in  Table  1  information  on  the  two  1'xl*  anomaly  fields  to  be  used 
in  our  combination  solution. 
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As  normal  figure  of  the  earth  a  geocentric,  rotating,  equipotential 
ellipsoid,  the  gravity  field  of  which  can  be  defined  by  four  parameters,  is 
formally  used.  The  parameters  can  be  the  equatorial  radius  a,  the  flattening 
f,  the  enclosed  mass  M  which  is  equal  to  the  earth's  mass,  and  the  rotational 
velocity  u  which  is  equal  to  the  earth’s  rotational  speed.  The  normal  gravity 
potential  generated  by  this  ellipsoid  can  be  expressed  in  terms  of  even-degree 
zonal  harmonics  as: 


U(r,  ♦  ,  A)  =  ^[l  +  JaDan>0Pan(sin*)]  +  |o2r2cos2i 
where 

D2n,o  fully  normalized  normal  potential  coefficients 
Pln  fully  normalized  Legendre  polynomials  of  degree  2n. 

The  D|  o  are  related  to  the  conventional  J|  coefficients  as  follows: 


(3.5) 


(3.6) 


21  +  1 


In  practice  the  J <  are  negligible  for  I  >  6.  From  Cook  'lOSi))  we  can  compute: 


j.  =  |  [fd  -  -  i  .(i  -  f  f  *  ii  f)j 

=  §i  m  -|f)[vf(i  -|f)  -  5.(1  -ffj) 


(3.7) 


(3.8) 


J.  =  2l  fJ(6f  "  5b,) 


(3.9) 


;2as(l-f) 


(3. 10) 


The  disturbing  potential  is  defined  as: 


T  =  W  -  U 


(3.11) 


Using  (3.2)  and  (3.5)  this  becomes: 


IfM  •  t  1  (all. 

T(r ,  ♦,  A)  =  T3  t  I  I  ?  C,*  Y|g(*.  A) 

r  l  =  a  m—  o  oi-o'r ’ 


(3.12) 


where  the  fully  normalized  disturbing  potential  coefficients  C|g  are  the  same 
as  the  A|g  except  in  the  case  of  the  even  degree  zonals  for  which  the  D|>0 
are  subtracted  from  the  A<(0.  In  the  usual  notation: 


(3.17) 
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(1 

+  i 

.  2cj2  e4rJsin224) 

T 

1  e4rssin22i 

Ih 

N 

y  4b 4M  J 

T  N 

ir  4b4 

where 

M  meridional  radius  of  curvature  of  the  ellipsoid 
N  pr  •  vertical  radius  of  curvature  of  the  ellipsoid 
b  polar  radius  of  the  ellipsoid 

Retaining  only  terms  of  0{e2)  the  boundary  conditions  (3.16)  and  (3.17) 
become  (Cruz,  1986,  (2.7),  (3.10)): 


Ag=-~-~T-  e2sin4cos$  — +  e2(3sin24  -  2)  — 
Jr  r  r3+  r 

agH  =  "  f?  "  \  T  +  e2(3sin24  -  2)  J 


(3.18) 

(3.19) 


If  terms  of  0(e2)  are  further  neglected  the  boundary  condition  transforms  to 
the  usual  spherical  form: 


9T  2 

Ag°  =  =  -  37  -  f  T 


(3.20) 


where  the  superscript  "o"  denotes  the  spherical  approximation  gravity 
anomaly. 


3.3  The  Gravity  Anomaly  in  Terms  of  Potential  Coefficients 


In  our  further  discussions  we  concentrate  on  the  use  of  the  gravity 
anomaly  Ag  defined  in  (3.14)  and  the  boundary  condition  (3.18),  as  opposed  to 
the  AgM  of  (3.15)  and  (3.19). 

We  first  wish  to  express  the  gravity  anomaly  Ag  of  (3.14)  in  terms  of  the 
disturbing  potential  coefficients  C occurring  in  (3.12).  This  is  conceptually 
done  by  substituting  (3.12)  into  (3.18).  The  result  is  formulated  in  Cruz 
(1986)  in  terms  of  ellipsoidal  corrections: 

Ag  =  7?  J  i  i  (i-i)(;)\c4s:  -  c?;h  -  C?;7)Y,£(i,  X)  (3.21) 

r  i  =  i  m—  o  ot  —  o  > 

where  the  C|^,'  (i  =  h,  y)  are  ellipsoidal  corrections  of  0(eJ)  computed  from 
the  C|°J  themselves  a3  follows: 

=  eJ(p  +  q|^C<«  +  rj,j,v*mC|  +  3 .  J ;  i  =  h,  y 


WM  m  t  i  ,  . 

*&°a  =  rr  ,1  I  I  (i-DCig  YiSSU,  x) 

«  1=2  m=o  <*=o 

with  inverse 

5  5*777-  s  JJ  ’**<♦•  x> d‘ 


13.28) 


(3.29' 


We  still  need  to  relate  the  tig°  needed  in  (3.29)  to  the  gravity  anomaly  that  is 
more  directly  available  in  practice,  namely,  the  gravity  anomaly  Ags  that 
follows  the  definition  (3.14)  and  refers  to  the  topographic  surface  of  the 
earth. 

To  relate  Ags  and  Agjj  Cruz  (1985)  has  suggested  that  Ags  be  downward 
continued  to  the  ellipsoid  as  a  first  step.  To  do  this  we  utilize  a  Taylor 
series  expansion: 


AgE  = 


_  IA&  h  _  i_  3.2Ag  h  2 
a r  2!  a r2 


( 3 . 30  'i 


where  h  is  the  height  of  the  point  (or  mean  elevation  of  a  block)  above  the 
ellipsoid.  In  writing  (3.30)  we  recognize  the  concern  with  analytical 
continuation  and  the  problem  of  computing  the  gradient  expressions.  Although 
local  gradients  can  be  computed  from  detailed  data,  we  need  gradients  for  the 
reduction  of  mean  anomalies  on  a  global  basis.  The  moBt  efficient  way  to  do 
this  is  to  differentiate  the  anomaly  expression  (3.21)  assuming  potential 
coefficients  are  known.  Neglecting  the  small  effects  of  and  in 

these  reductions  we  have: 


IAS  - 

3r 


ifM  ®  4  i 

•  3  i  l  l  (i-l)(<+2)  f  X) 

r  1-2  m=o  tx-o  'T  > 

“  1  i  i  (<-l)(i+2)(i+3)(2)iCjS<  Y|«fi,  X) 

r*  £-2  m - o  c<  =  o 


(3. 31' 


(3.32) 


The  root  mean  square  value  of  the  first  derivative  reduction  in  (3.30)  is  *10 
ragal  (max  =  19  mgals)  and  for  the  second  derivative  reduction  it  is  *0.01  mgal 
(max  c  1.0  rngal)  (Rapp,  1984).  These  computations  were  carried  out  with  the 
Rapp  (1981)  potential  coefficient  field  to  degree  180.  From  recent  calculations 
using  l*xl*  mean  gradient  computation  and  an  unpublished  expansion  to 
degree  300,  the  RMS  value  of  the  first  derivative  effect  is  *0.8  mgal  (max  =  30 
mgals)  and  for  the  second  derivative  it  is  *0.04  mgal  (max  =  2.1  mgals).  In 
practice  we  compute  the  gradient  correction  terms  and  store  them  for 
application  to  various  sets  of  surface  gravity  anomalies. 

Knowing  the  gravity  anomaly  Agf  on  the  ellipsoid  we  can  express  the 
spherical  approximation  gravity  anomaly  Ag£  on  the  equatorial  sphere  us  (Cruz, 
1986,  (4.8)): 


=  AgE  +  'b‘‘,7*£i  +  »2  +  '3* 


'3.33' 
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and 


Ci"’J  =  sTTTT  F(1-5)(<-2)(<-i)Ki-.,mwimc?i:,(1, 
+  (i-3)(l)(<  +  l)L<-2tn)ulmC?i“(m 
+  (4-l)(4+2)(4+3)M<mC?;° 

+  (4+l)(4+4)(4+5)Ni+2)mvimCi;°im 

+  ( *+3)  (*+6)  (<+7)Pj+4 f  (Bx|mCj4-4,m] 

with 

*■  ^im^l+2  ,fn 

^*lm  ~  +  P  & +■  2fm) 

~  +  2,m  *  —  2,m7im 

"  7 1  rr\  fi  t  ~  2  ,m  +  Pitn^ 

~  7i~ 2 , mTlm 

(*-m+l)(l-m+2) 

°'<m  "  (24+l)(24+3) 

(24*~2ma+24-l) 

~  (2<+3)(2/-l) 

( 4  +m)  (  4 +m~l ) 

ytm  -  (2I  +  1H24-1) 


(24~7)(4+ro-3)(4+a-2)(4+m+l)(4+i) 
w<m  *  y  (24  + 1 )  ( 4-m~3)  ( 4-m-2) ( 4-n-l )  ( 4-n) 

f  (24+9)(4-m+-l)(4-m+2)(4-in+3)(4-B>+-4) 
X<m  "  /  ( 24+ 1 )  ( 4+m+l )  ( 4+n»+2 )  (  4+B+-3  )  ( 4  +  m+-4  ) 


! 


(3.39) 


v3.40' 


(3.41  > 


3.  12' 


This  completes  the  discussion  on  how  to  express  the  disturbing  potential 
coefficients  C|{*  in  terms  of  the  surface  gravity  anomaly.  The  two  steps 
involved  are  (3.30)  followed  by  an  implementation  of  (3.34)  via  (3.35)  or  (3.36). 
the  formula  for  the  Cj£>3  in  terms  of  Cj£'°  is  given  in  Cruz  (1986,  (6.20)). 


3.5 _ Physical  Realization  of  the  Gravity  Anomaly 

A  final  point  should  be  made  concerning  the  physical  realization  of  the 
theoretical  gravity  anomaly  defined  in  (3.14).  Rapp  (1984)  discusses  the 
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The  mathematical  model  is  written  as 

F  =  F(Lia,  L 

«■)  =  o 

(4.1) 

which  is  linearized  to  yield  the  observation  equation: 

B|V<  ♦  BXVX 

+  W  =  0 

(4.2) 

where 

„  a  f 

„  a  f 

Bl  '  3Li: 

*  =  *LX 

(4.3) 

If  we  designate 

Pi  and  Px  as  the  weight  matrices  for 

the  observations 

parameters  respectively,  the  weighted  least  squares  condition  for  the  solution 
is: 


ViTP|V|  +  VjPxVx  =  a  minimum  (4.4) 

The  solution  for  Vx  is: 

Vx  =  -(BjM~'Bx  +  Px)BjM~‘W  (4.5) 

where 

M  =  BfP| 1 B| 1  (4.6) 

The  observation  residuals  are: 

V,  =  -P*‘B*TM  l(BxVx  +  W)  (4.7) 

The  error  variance -covariance  matrix  for  the  solution  vector  is: 

£x  =  mo (BjM“lBx  +  Px)-‘  (4.8) 

where  m^  is  the  variance  of  unit  weight. 


The  mathematical  structure  is  formed  by  assuming  that  the  potential 
coefficient  estimate  from  satellite  derived  procedures  should  be  the  same  as 
that  obtained  from  the  global  gravity  data  set  after  appropriate  corrections 
(e.g.  ellipticity  and  downward  continuations)  have  been  made.  We  write: 

F  =  Lxa  -  Lx a  =  0  4.9) 

r 

If  Lxo  is  the  approximate  potential  coefficient  set  and  Lxc°  is  the  coefficient 
set  based  on  the  observed  (and  corrected)  anomalies,  the  misclosure  vector  is: 

W  =  Lxo  -  I,  c  (4.  10) 

*0 

The  Lxc  values  are  computed  from  (see  eq.  3.34) 

c,i'°  "  1^,-1)  II  *«>**«<♦•  x>  **  l4-11' 

<y 
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The  corresponding  adjusted  potential  coefficients,  considering  the  ellipsoidal 
corrections,  would  be: 


i 


Ci, 


Lx*  *  6Ct«’°  + 


The  adjusted  anomalies  would  be  found  from: 


Lfa  =  L |°  +  Vi 


4.21 


(4.22) 


The  corresponding  surface  anomalies  would  be  found  from  (4.13) 


L«a 


+ 


3Ag 

3r 


h  + 


i_  Has 

2!  3r2 


h* 


'4.23' 


This  completes  the  discussion  of  the  adjustment  model  to  be  used  in  this 
combination  solution.  It  is  similar  to  the  procedure  used  in  the  computation  of 
the  OSU78  combination  solution  except  now  we  reduce  surface  anomalies  to  the 
ellipsoid  and  apply  ellipsoidal  corrections. 


The  normal  equation  matrix  of  the  solution  is  given  by  BiPj'Bf1.  The  size 
of  the  matrix  is  dependent  on  the  number  of  a  priori  defined  potential 
coefficients.  For  a  solution  complete  to  degree  20  the  number  of  coefficients 
is  441.  In  the  1978  solutions  we  were  able  to  complete  a  solution  to  t  -  12. 
In  our  1981  solution  we  wanted  to  go  to  a  much  higher  degree  with  more 
coefficients.  This  could  not  be  done  because  of  computer  limitations  so  we 
made  certain  assumptions  that  led  to  the  normal  equation  matrix  being 
diagonal.  The  result  was  an  approximate  combination  solution  that  was  carried 
out  with  small  computer  demands. 


For  the  solution  to  be  described  in  this  report  we  have  64800  l*xl* 
"observations"  and  582  unknowns.  The  estimated  computer  time  for  a  solution 
on  the  IBM  4031D  at  OSU  was  4.5  hours.  As  several  solutions  were 
contemplated  the  OSU  computer  requirements  were  still  considered  excessive. 
At  this  point  we  started  tests  on  a  CRAY  XMP2/4  supercomputer  at  Boeing 
Computer  Services.  Funding  for  this  computer  work  was  from  National  Science 
Foundation  Grant  EAR-8420862.  Modification  of  our  existing  combination 
program  by  one  of  us  (JC)  took  advantage  of  the  vectorization  ability  of  the 
CRAY  machine.  With  this  change,  a  combination  solution  took  6.5  minutes. 
This  meant  that  a  number  of  different  solutions  could  be  carried  out  to  learn 
what  the  best  solutions  could  be. 


5.  The  Batimation  of  the  High  Degree  Potential  Coefficients 

The  adjustment  described  in  the  previous  section  yields  a  set  of  adjusted 
potential  coefficients  (see  eq  4.20)  and  a  set  of  adjusted  anomalies  on  the 
ellipsoid  (see  eq.  4.22).  These  anomalies  can  then  be  used  to  estimate  Cf£>° 
using  the  following  (see  eq.  3.34) 


CiSS>° 


7~  II  Agf  YjS${ ♦ ,  A)  d<r 


;5.  i 


where  N  =  18O*/0*.  The  evaluation  of  (5.2)  using  pre-computcd  associated 
Legendre  function  integrals,  and  a  Fourier  analysis  of  the  mean  anomalies 
along  latitude  bands,  and  with  the  q|  values  defined  in  (5.5)  is  carried  out 
through  the  HARMIN  subroutine  in  Colombo  (ibid)  as  implemented  in  OSU  (RHR) 
library  program  F419.  This  procedure  is  one  of  the  ways  in  which  the  high 
degree  expansions  have  been  computed  from  the  adjusted  anomalies. 

5.2  Potential  Coefficients  Through  Ontimal  Estimation 

A  disadvantage  of  the  procedure  in  the  previous  section  is  that  it  does 
not  take  into  account  data  noise  nor  does  it  provide  error  estimates  of  the 
estimated  coefficients.  To  consider  this  factors  we  have  used  the  optimal 
estimation  procedure  developed  by  Colombo  (1981)  and  implemented  for  l*xl* 
data  by  Hajela  (1984).  We  discuss  only  the  general  principles  of  this  least 
squares  collocation  solution  to  (5.1).  We  start  with  the  global  mean  gravity 
anomaly  vector,  dg,  expressed  as  the  sum  of  a  signal  vector,  z,  and  a  noise 
vector,  n: 

dg  -  z  +  n  (5.6) 

Let  the  potential  coefficient  vector  (Cfj**0]  be  defined  as  c.  Let  F  be  a 
linear  operator  that  will  estimate  c  (i.e.  c)  from  dg: 

c  -  F(z  +  n)  (5.7) 

The  error  in  this  estimation  is  e: 

§=c-c-c-  F(z  +  n)  -  (c  -  Fz)  -  (Fn)  '5.B' 

We  define  the  sampling  error  as  es  and  the  propagated  noise  error  to  be  e„ . 
We  have: 


es  -  c  Fz  ;  en  -  Fn  .5.9 

The  error  covariance  matrix  of  e  is: 

et  r  f  En  =  M(es  e^}  +  M{en  ej}  ‘5.10: 

where  M  is  an  averaging  operator.  Substituting  (5.9)  into  (5.10)  we  find: 

E,  =  C  -2CczFt  +  F(CZ7  +  D)FT  (5.11 

where 

C  =  M{c  cM,  Ccz  --  Mfc  z T } 

5.12' 

=  M{z  z1},  D  =  M{n  nT) 

C  represents  the  covariance  matrix  of  the  potential  coefficients;  Ccl  is  the 
cross  covariance  between  the  potential  coefficients  and  the  given  mean 


The  error  variance  of  the  potential  coefficients  is  the  same  for  C1  and  C2 . 
It  is  given  by: 


a  c  t  m  ~  ti\ 


1 


72(t-D2 


_£t_ 

2<  +  I 


AX 2 

1-cosmAX 


2N  xii 

i  =o 


if  m  =  0 
if  m  *  0 


(5. 18) 


where  c*  are  the  anomaly  degree  variances  implied  by  an  a  priori  potential 
coefficient  model. 

The  calculation  of  the  covariance  functions  and  the  potential  coefficient 
errors  of  the  optimum  estimation  solution  requires  the  definition  of  c|. 
Adopted  here  was  a  model  used  by  Colombo  (ibid)  which  consists  of  the 
anomaly  degree  variances  implied  by  the  OSU78  solution  from  degree  2  to  100 
and  beyond  100  by  a  two  component  anomaly  degree  variance  model  described 
by  Rapp  (1979).  Specifically  this  model  (with  the  c|  values  interpreted  to  be 
at  a  sphere  of  radius  R  (=  6371  km)  is: 


ci(R) 


+ 


1 

i+B  1-2 


(s2)*  +  2]  nagals 


•  5.  19' 


where 

o i,  =  3.4050,  cxj  =  140.03 
s,  =  0.998006  s  2  =  0.914232 
A  =  1,  B  =  2 


For  actual  implementation  we  require  the  c|  values  on  a  sphere  of  radius 
a  so  we  used: 


c*(a) 


1+j 


(5.20' 


A  plot  of  this  function  will  be  given  in  a  later  section  dealing  with  results. 

Our  intention  in  this  study  will  be  to  implement  both  the  HARMIN  type 
solution  and  the  optimal  estimation  solution.  Various  solutions  will  be  made 
and  the  results  compared. 


6.  Preliminary  Investigations 

Before  we  discuss  our  final  solutions  it  is  appropriate  to  describe  a  set  of 
studies  that  led  to  the  procedures  used  in  the  final  solution.  These  studies 
were  carried  out  with  preliminary  sets  of  global  gravity  models  and  a  priori 
potential  coefficient  sets. 
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Combination  Solutions 


R.  Cumulat  ively: 


,6.1; 


(2)  The  root  mean  square  anomaly  difference  by  degree  and  for  the  whole 
coefficient  set: 


A.  By  degree 

-  l?2(t  I)1  l  1  acJ;*P 

1  m~0  ct-o  J 

6.3' 

B.  Cumulatively: 

f  Nma  x  ] \ 

<5t?  =  1  I  <5g,J 

1  1  -  2  1 

6.4 ; 

(3)  The  percentage  difference  by  degree  and  cumulatively: 
A.  By  degree: 


Pi 


I  I  AC, 

■  jP  -  Q 


I  I  cj;* 

n-o  al- 0 


100 


B.  Average  Percentage  Difference 


Nroa  * 

X  Pt 

p  -  L.7-1 _ 

'  Nma „ - 1 


6. 5 


6.  H' 


The  valueH  of  <5N  and  of  6g  represent  the  global  mean  difference  of  the 
two  potential  coefficient  functions  over  the  sphere.  The  values  of  6N,  and 
fig,  represent  the  same  information  by  degree.  The  percentage  difference 
provides  information  that  takes  into  account  the  relative  magnitude  of  the 
coefficient  differences  to  the  coefficient  magnitudes. 


Another  quantity  that  can  be  used  for  the  comparisons  of 
sets  is  the  degree  correlation  coefficients  (pi)  and  the  overall 
coefficient  p.  The  degree  correlation  between  the  coefficient 
would  be: 


Pi 


l  (Ct 

m-  Q 


m  j  ^  tc 


S|m  ) 


[  I  i  ]%  (  1  i  c If  ]■ 

Lm=o  o(=o  1  1  lm-o  «=o  j  J 


two  coefficient 
correlation 
sets  i  and  j 


(6.7; 


The  overall  correlation  would  be: 


p 


"man  t 

X  l  Wtmct_ 

l=j  m— O  1  i 


["f  r  i  ci:“  ]'  ["!"  1  i  cji-  ] 

'■t~  2  m-o  or-  o  '  1  '■1-2  m-o  0 1-0  J  ‘ 


(6.81 


variances  implied  by  a  uniform  uncorrelated  error  of  *10  mgals  in  l*xl* 
blocks.  From  Rapp  (1981),  we  have  that  the  propagated  standard  deviation  of 
a  fully  normalized  potential  coefficient  of  degree  t  due  to  a  uniform, 
uncorrelated  error  of  m(dg)  in  a  block  size  of  0  (radians)  can  be  approximated 
as: 

m(C,  S)  =  — 1E1A£2A  (6.9) 

The  10  mgals  used  for  m(Ag)  in  (6.9)  is  the  root  mean  square  (RMS)  of  the 
8-15  mgal  standard  deviations  used  in  the  combination  solution.  The  similarity 
between  columns  (3)  and  (5)  is  evident. 

To  examine  whether  the  assumption  of  8-15  mgal  uncorrelated  errors  was 
reasonable  we  computed  a  set  of  potential  coefficients  from  only  the 
unadjusted  l*xl*  anomalies.  These  values  were  then  differenced  from  the 
unadjusted  GEML2'  coefficients.  The  difference  anomaly  degree  variances 
were  then  computed  and  shown  in  Table  5,  column  (2).  It  is  seen  that  the 
data  errors  shown  in  columns  (1)  and  (3)  cannot  account  for  the  differences 
shown  in  column  (2).  This  implies  that  the  8-15  mgal  uncorrelated  noise 
assumption  is  very  optimistic. 

We  decided  to  examine  the  use  of  an  RMS  standard  deviation  of  *25  mgals 
for  the  l*xl*  anomalies.  The  assumption  of  uniform  uncorrelated  errors  leads 
to  error  anomaly  degree  variances  that  are  Bhown  in  Table  5,  column  (4). 
These  variances  are  now  (25/10) 2  =  6.25  times  larger  than  those  of  the  *10 
mgal  accuracy,  and  appear  to  be  much  more  reasonable.  The  *25  mgal  RMS 
accuracy  could  be  implemented  in  our  combination  solution  by  multiplying  the 
8-15  mgal  accuracies  by  2.5  then  restricting  the  range  to  20  «  m  ‘  38  mgals. 
We  then  performed  a  combination  solution  denoted  by  CS(20,  38),  yielding 
error  anomaly  degree  variances  for  the  adjusted  coefficients  as  shown  in 
Table  5,  column  (6).  These  values  are  now  larger,  with  cumulative  variance  to 
degree  20  of  1.88  mgal2  compared  with  0.48  mgal2  from  column  (5). 

It  is  of  interest  to  compare  the  coefficients  of  the  CS(8,  15)  and  CS(20,  38) 
solutions  in  terms  of  undulation  and  anomaly  differences  and  percentage 
differences.  Using  the  equations  of  section  6.2  the  difference  between  the  two 
solutions  are  shown  in  Table  6.  As  would  be  expected  the  greatest  change 
takes  place  at  the  lower  degrees  because  it  is  here  that  the  dual  estimate  of 
the  potential  coefficients  are  obtained  from  the  adjustment  of  satellite  and 
terrestrial  data. 

The  CS(8,  15)  and  CS(20,  38)  types  of  solution  were  also  compared  in 
terms  of  orbit  fits  and  fits  to  Doppler  derived  undulations.  The  orbit  fit 
results  are  shown  in  Table  7,  with  more  specific  details  of  orbit  tests  given  in 
Section  8.5.  It  is  seen  that  the  CS(20,  38)  solution  performs  substantially 
better  than  the  CS(8,  15).  This  seems  natural  as  the  higher  anomaly  standard 
deviations  gives  the  a  priori  potential  coefficients  greater  influence  on  the 
solution. 


Table  8.  Comparison  of  Doppler  Derived  Undulations  with  Preliminary 
Geopotential  Models  Complete  to  Degree  180. 


Mean 

Std  Dev. 

Number  of 

Difference 

Difference 

Stations 

Area 

CS(8, 15) 

0.32  m 

*1.64  m 

1749 

Global1 

CS(20,38) 

0. 19 

1.68 

1738 

CS(8,  15) 

0.38 

1.55 

691 

N.  America1 

CS(20,  38) 

0.21 

1.54 

687 

CS(8, 15) 

-0.24 

1.42 

172 

Europe2 

CS(20,  38) 

-0.48 

1.38 

172 

CS(8,  15) 1 

-1.13 

1.60 

114 

Australia2 

CS(20,  38) 

-1.06 

1.69 

114 

1  stations  with  residuals  greater  than  4  meters  in  absolute  value  were 
rejected 

2  A  fixed  set  of  stations  was  used 

J  The  combination  solution  used  an  anomaly  standard  deviation  range  of 
8  <  m  <  15  mgals  over  the  land  of  Australia  only,  with  20  <  m  *  38  mgals 
used  over  the  rest  of  the  world. 


difference  standard  deviation  of  *1.60  m  at  114  Australian  Doppler  stations 
with  our  final  solution  to  be  discussed  later,  compared  with  *1.69  when  the 
CS(20,  38)  solution  was  used. 

The  problems  discussed  in  this  section  arise  from  the  assumption  of 
uncorrelated  noise.  In  fact  the  anomaly  errors  are  correlated.  Such 
correlation  and  its  role  in  the  computation  of  anomaly  degree  variance 
accuracies  is  discussed  by  Weber  and  Wenzel  (1982).  At  this  point  we  have  no 
way  of  treating  anomaly  error  correlation  in  our  solutions.  Instead  we  have 
increased  our  anomaly  error  estimates.  Such  a  procedure  is  not  desirable  and 
in  fact  gives  unreasonably  high  errors  at  degrees  above  30.  To  compensate 
for  this,  our  accuracy  estimate  for  unadjusted  coefficients  will  be  based  on 
the  *10  mgal  uncorrelated  noise  assumption. 


6.4  Preliminary  High  Degree  Expansions 

We  next  studied  the  manner  by  which  the  high  degree  expansions  from  the 
adjusted  anomalies  of  the  combination  solution  should  be  carried  out.  The 
main  comparisons  to  be  discussed  here  relate  to  a  HARMIN  type  solution  (see 
equation  5.2)  and  the  optimal  estimation  solution  (see  equation  5.15).  The  first 


implied  by  the  0E{8,  15)  and  0E1  test  solutions  as  well  as  the  a  prion  model 
used  in  the  optimal  estimation  procedure.  The  smoothing  effect  of  the  anomaly 
noise  is  evident  from  this  plot. 

The  optimal  estimation  procedure  yielded  error  estimates  for  the  computed 
coefficients.  For  the  0E(8,  15)  solution  these  estimates  were  directly  taken, 
giving  error  anomaly  degree  variances  as  also  plotted  in  Figure  7.  For  the 
0E1  solution  the  error  estimates  from  the  optimal  estimation  mainly  correspond 
only  to  the  sampling  error,  with  a  small  contribution  from  the  1  mgal  noise. 
To  complete  these  estimates  we  quadratically  added  the  propagated  error 
implied  by  the  assumption  of  *10  mgal  uniform,  uncorrelated  l*xl*  noise. 
Specifically,  we  computed: 


mfC,  S)  =  v  mJ  f  SE ) 


16.11 


where  m(SE)  is  the  sampling  error  and  m(PE)  is  the  propagated  error.  The 
m(PE)  is  given  by  (6.9).  A  similar,  but  not  as  rigorous  procedure  was  used  in 
the  development  of  the  OSU81  field.  The  total  error  estimate  for  the  0E1 
solution  gave  error  anomaly  degree  variances  that  are  also  shown  in  Figure  7. 


The  error  estimates  for  the  OE(8,  15)  solution  are  smaller  than  the  total 
error  estimates  of  the  OE1  solution.  This  is  consistent  with  the  concept  of 
filtering  the  anomaly  errors  in  the  OE(8,  15)  solution.  However,  this  filtering 
involved  a  smoothing  of  the  spectrum  which  could  substantially  affect  valid 
information.  In  order  to  test  the  two  solutions  we  compared  then  in  terms  of 
their  fits  to  Doppler  derived  undulations  as  shown  in  Table  10.  The  OH) 
solution  performed  slightly  better  in  these  comparisons.  At  this  point  we 
decided  to  use  for  our  work  the  0E1  type  of  solution,  with  a  simple 
propagation  of  the  anomaly  errors  as  shown  in  (6.11). 


Since  the  optimal  estimation  solution  is  a  complicated  process  it  is  of 
interest  to  compare  the  results  of  HARMIN  with  OE1.  These  comparisons  are 
shown  in  Table:  11.  The  conclusion  from  Table  11  is  that  HARMIN’  gives 
coefficients  that  agree  well  with  the  coefficients  found  from  the  optimal 
estimation.  The  disadvantage  of  the  HARMIN  solution  is  that  no  error 
estimates  are  provided.  Such  values  are  found  when  the  optimal  solution  is 
carried  out. 
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Table  10.  Comparison  of  0E(8,  15}  and  OE1  in  Terms  of  Fits  to  Doppler 

Derived  Undulations.  Maximum  Degree  is  180. 


Mean 

Std  Dev. 

Number  of 

Model 

Difference 

Difference 

Stationst 

Are  u 

0E(8, 15) 

0.22  m 

*1.68  m 

1731 

Gl ohal 

0E1 

0.  19 

1.68 

1738 

0E(H,  15) 

0.25 

1.56 

683 

N.  America 

0E1 

0.21 

1.54 

687 

t  Stations  with  residuals  greater  than  4  meters  in  absolute  value  were 
rejected 


Table  11.  Comparison  of  Potential  Coefficient  Solutions  From  HARMIN  arid  from 
Optimal  Estimation  (1  mgal  standard  error). 


t 

6N  (cm) 

<5g  (mgals) 

P(\) 

2 

0.0 

0.00 

0.0 

10 

0.0 

0.00 

0.0 

20 

0.0 

0.00 

0.0 

30 

0.2 

0.01 

0.5 

50 

0.2 

0.02 

1.1 

75 

0.8 

0.09 

6.2 

100 

1.0 

0.  15 

10.6 

120 

1.0 

0.  18 

13.2 

150 

0.6 

0.  14 

13.0 

180 

0.8 

0.22 

22.3 

to  1  HO 

8.9 

1.53 

7.6 

6.5  A  Study  of  the  Undulation  Residual  Correlation  with  Elevation 

Tacherning  (1985,  private  communication)  pointed  out  that  the  differences 
in  the  Doppler  and  gravity  model  undulations  were  correlated  with  elevation 
for  the  0SU81  solution  but  not  for  the  GPM2  field.  In  our  research  we  have 
computed  this  correlation  for  several  models  for  various  Doppler  data  sets  by 
assuming  there  was  a  linear  correlation  with  elevation.  To  do  this  the 
available  Doppler  stations  were  grouped  into  100-meter  elevation  intervals. 
The  undulation  residuals  AN  llirig  within  an  interval  were  then  averaged  to 
get  AN.  The  mean  value  AN  was  assumed  to  refer  to  the  mean  elevation  H  of 
the  stations  involved.  A  simple  least  squares  fit  was  then  carried  out  to 
compute  the  slope  s  and  bias  b  in  the  observation  equation: 
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Table  12.  Undulation  Residual  Correlation  with  Elevation  Using  the  Global 
Station  Set.  Maximum  Degree  is  180. 


Solution 

Slope 

(m/km) 

Std.  Dev. 
(m/km) 

Number  of 
Stat lonst 

0SU81 

0.73 

*0.16 

1721 

GPM2 

0.00 

0.10 

1735 

CS(20 ,  38,  0E1) 

0.39 

0.  13 

1738 

GPM2(n  =  7  to  14) 

+  CS(20,  38,  OF  1 ) 

0.02 

0.  13 

1734 

CS( 8,  15,  OKI) 

'1.42 

0.  13 

1749 

t  stations  with  residuals  greater  than  4  meters  in  absolute  value  w.-t- 
rejected 


The  slopes  of  our  new  solutions  are  less  than  in  OSt'Rl  ('0.4  rn/km 
compared  t.o  '0.7  m/km),  although  they  may  still  be  considered  sign-ficant. 
Further  teats  are  needed  to  better  understand  the  apparent  slope  pr  and 

its  implications. 


Table  1.1.  Undulation  Residual  Correlation  with  Elevation  1  sing  th<  Ikdi.il 
Station  Set.  Maximum  Degree  is  30. 


Sol ut ion 

S  1  ope 

m  km ' 

Std.  Dev. 

m  km 

Number  of  , 
St  at  i onst 

1 

OS  UR  1 

0.93 

*0.  28 

1694  1 

GPM2 

0.24 

0.17 

1725  ; 

GFML2 

0.57 

0.21 

16  19  j 

GRIM3I.1 

0.79 

0.19 

It. 92  : 

OS ' 20 ,  38,  OF l 

0.74 

0.  17 

it:u  j 

_ i 

+  Stations  with  residuals  greater  than  5  meters  in  absolute  value  wen- 
reject  ed. 


t 

( 
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In  addition  we  found,  again  through  Doppler  undulation  comparisons,  some  1 

geophysical  anomalies  to  be  quite  realistic  in  the  sense  that  they  yield  , 

improved  undulation  values  at  selected  sites.  Because  of  this  we  selected  (see 
Section  6.6)  63  geophysical  anomalies  to  be  included  in  the  data  set  for  which 
almost  all  geophysical  anomalies  were  excluded.  The  formation  of  this  global 
data  set  was  discussed  in  Section  2.4. 

Tests  with  these  fields  indicated  that  our  Doppler  undulation  comparisons 
in  Australia  were  not  as  good  as  found  in  other  areas,  primarily  North  America 
and  Europe.  In  order  to  obtain  a  better  Doppler  fit  in  the  Australia  region 
we  gave  the  terrestrial  anomaly  data  higher  weight  by  specifying  in  this 
region  the  anomaly  standard  deviation  range  to  be  8  to  15  mgals  instead  of 
the  20  to  38  mgal  used  everywhere  else  (see  Section  6.3).  In  addition  we 
changed  five  anomalies  from  one  data  source  to  another  source  on  the  basis  of 
individual  Doppler  undulation  comparisons.  This  led  to  undulation  differences 
decreasing  from  *1.69  m  to  *1.60  m.  More  details  of  undulation  comparisons 
are  given  in  Section  8. 

With  the  above  as  background  we  carried  out  the  two  main  solutions  of 
this  report,  OSl'86C  and  OSU86P  The  C  solution  essentially  excludes 
geophysical  anomalies  while  the  D  solution  includes  them.  The  combination 
solutions  were  performed  on  a  CRAY  XMP  2/4  machine  using  the  procedures 
outlined  in  our  previous  sections.  The  ellipsoidal  corrections  were  applied  to 
the  adjusted  coefficients  (see  eq.  4.21)  to  obtain  our  final  "adjusted" 
coefficients.  The  adjusted  anomalies  were  developed  into  potential  coefficients 
using  the  optimal  estimation  procedure  with  a  uniform  anomaly  noise  of  *1 
mgal  and  the  a  priori  anomaly  degree  variance  model  defined  by  equations  5.19 
and  5.20.  Th»  expansion  was  carried  out  to  degree  250  which  involves  the 
determination  of  63,001  coefficients.  The  coefficients  were  merged  with  the 
adjusted  coefficients  (see  remark  below  eq.  5.1),  ellipsoiaal  corrections 
computed  from  the  merger  (eqs.  3.37  and  3.39),  and  these  corrections  applied 
(eq.  3.36)  to  obtain  the  final  high  degree  coefficient  set.  Although  the 
adjusted  eoeff  icients  and  the  coefficients  computed  from  the  adjusted 
anomalies  should  be  the  same,  we  gave  preference  to  the  directly  adjusted 
potential  coefficients  because  such  coefficients  carried  a  standard  deviation 
determined  through  the  adjustment. 

The  standard  deviations  for  the  potential  coefficients  not  estimated  in  the 
adjustment  were  obtained  considering  the  sampling  error  (SE),  found  from  the 
optimal  estimation  solution  (see  eq.  5.18)  and  a  propagated  error  (PE)  implied 
by  the  assumption  of  *10  mgal  uncorrelated  noise  for  the  l’xl*  anomalies. 

This  error  is  computed  from: 


n (  PE) 


an  Agj B 


(7.1) 


2y  '  I  1  rr 

where  0  is  the  block  size  (I*  in  radians)  and  m(Ag)  is  the  anomaly  standard 
deviation  in  units  of  y.  The  coefficient  error,  for  either  C  or  S  was  then: 


m '  C  ,  S  )  ■  y~ m 2  (  SE  ~  m 2  (  PE ) 


(7.2) 
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8.  Comparison  of  Results 

The  purpose  of  this  section  is  to  compare  the  new  models  between 
themselves  and  with  two  other  high  degree  models.  The  models  are  the  0SU81 
(Rapp,  1981)  and  the  GPM2  (Wenzel,  1985)  model.  The  comparisons,  carried  out 
in  several  different  ways  are  described  in  the  following  section. 

8.1  Anomaly  and  Undulation  Comparison 

The  global  differences  between  the  various  solutions  are  shown  in  Table 
15  for  undulations  and  Table  16  for  anomalies.  The  differences  are  given  for 
the  maximum  degrees  of  30  and  180. 


Table  15.  Potential  Coefficient  Differences  in  Terms  of  Geoid  Undulations 
(meters) 


OSH  Dec:81 

to  t  =  30 
to  i  =  180 


0SU86D 

0SUB6C 

GPM2 

*  1.04 

*  1.32 

*  1.12 

1.20 

1.51 

1.31 

*  1.08 

1.21 

1.25 

1.43 

*  0.56 

0.73 

0SU86C 
to  i 
to  i 


Table  16.  Potential  Coefficient  Differences  in  Terms  of  Gravity  Anomalies 
(mgals) 


0SU86D 

0SU86C 

GPM2 

OSl'  Dec81 

to  t  -- 

30 

*  2.5 

*  3.0 

*  2.6 

to  t 

180 

7.3 

8.9 

8.4 

GPM 2 

to  *  = 

30 

*  2.6 

3.0 

to  i  - 

180 

8.5 

9.9 

0SU86C 

H 

O 

+-> 

30 

*  1.5 

to  t  - 

180 

4.9 

Undulation  Difference  Map  -  OSU86D  Minus  GPM2  (Contour  Interval 


We  next  turn  to  a  comparison  of  the  geoid  undulations  derived  from  (he 
various  geopotential  models  with  undulations  derived  from  Doppler  derived 
positions.  If  h  is  the  ellipsoidal  height  of  a  station  derived  through  Doppler 
satellite  positioning  techniques,  after  conversion  to  a  geocentric  system 
properly  scaled  (Rapp,  1983),  the  "Doppler  undulation"  is: 

N0  -  h0  -  H  8.1 

where  H  is  the  orthometric  height  of  the  point.  We  compute  N  from  Brim's 
equation  where  the  disturbing  potential  is  given  by  equation  (3.12)  with  r  the 
geocentric  distance  to  the  point  in  question  projected  to  the  ellipsoid.  In  all 
our  comparisons  we  used  the  parameters  to  convert  the  Doppler  system  to  a 
geocentric  system  with  correct  scale  derived  by  Boucher  and  Altarnimi  (19Hf>). 
The  translation  and  scale  parameters  to  go  from  the  Doppler  system  to  a 
geocentric  system  were  used  as  follows: 

Ax  -  10.  6  cm 

Ay  -  69.7  cm 
A 7  -  490 . 1  cm 
As  -  0.604  ppm 

The  comparisons  made  were  with  two  primary  sets  of  Doppler  stations. 
One  contained  approximately  800  stations  and  was  received  from  Tscherning. 
This  set  was  originally  obtained  from  the  National  Geodetic  Survey  and 
contains  stations  only  in  North  America.  The  second  set  contained 
approximately  2000  stations  distributed  globally.  The  coordinates  wen- 
determined  in  the  years  from  1971  to  1985.  The  number  of  passes  could  range 
from  a  minimum  of  25  to  a  maximum  of  592.  In  our  analysis  we  made  no 
correction  for  sun  spot  effects  as  suggested  by  Tscherning  and  Goad  (1985). 
The  overall  effect  of  neglecting  this  for  our  large  data  sets  is  expected  to  be 
only  a  few  cm.  In  making  our  comparison  we  assumed  the  geoid  undulations 
referred  to  an  ellipsoid  with  equatorial  radius  of  6378136  m  and  a  flattening  of 
1/298.257.  Comparisons  for  the  Tscherning  data  set  are  given  in  Table  19. 
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8.4  Undulation  Residual  Correlation  with  Elevation 


Tscherning  (1985,  private  communication)  pointed  out  that  the  differences 
in  the  Doppler  and  gravity  model  undulations  were  correlated  with  elevation 
for  the  0SU81  solution  but  not  for  the  GPM2  field.  In  our  research  we  have 
computed  this  correlation  for  several  models  for  various  Doppler  data  sets  by 
assuming  there  was  a  linear  correlation  with  elevation.  Results  for  the  global 
data  set  and  the  stations  in  Europe  are  shown  in  Table  23  and  24. 


Table  23.  Undulation  Residual  Correlation  Using  the  Global  Station  Set 


Solut ion 

Slope 

(«/ka) 

Std  Dev. 

(m/km) 

No .  of 
Stations 

0SF81 

0.73 

*  0.16 

1721 

GPM2 

0.00 

0.10 

1735 

0SL186C 

0.41 

0. 12 

1741 

0SU86D 

0.43 

0. 12 

1743 

Table  24.  Undulation  Residual  Correlation  Using  the  European  Station  Set 


Solution 

Slope 
; m/km) 

Std  Dev. 
(m/km) 

No.  of 
Stat ions 

0SU81 

1.53 

*  0.81 

172 

GPM2 

1.87 

0.63 

172 

0SU86C 

0.92 

0.65 

172 

0SU86D 

1.14 

0.69 

172 

From  these  results  we  see  that  the  GPM2  model  yields  residuals  with  no 
elevation  correlation  for  the  global  data  set  but  significant  correlation  for  the 
European  stations.  Tests  described  in  Section  6.5  indicate  that  the  correlation 
in  the  non-GPM2  solution  arises  from  the  gravity  field  information  in  degrees 
7  thru  14.  Further  tests  are  needed  to  better  understand  the  apparent 
correlation  and  its  significance. 


8.5  Correlation  with  Topography 


Under  certain  assumptions  the  topography  and  its  isostatic  compensation 
can  be  considered  to  generate  a  portion  of  the  earth's  gravity  field.  It  is 
appropriate  to  consider  the  correlation  of  the  spherical  harmonic 
representation  of  the  topographic/isostatic  potential  and  of  the  gravitational 
potential  represented  by  various  models.  Such  correlation,  by  degree,  can  be 
computed  from  eouation  (6.7).  The  average  correlation  for  several  different 
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Table  26.  Root  Mean  Square  Residual  From  Satellite  Orbit  Fits  Using  Selected 
Geopotential  Models 


Potential 

Model 

Satel 1 ite 

Standard* 

86C 

86D 

GPM2 

Starlette 

1  47  cm 

143 

179 

737  cm 

Seasat 

*  0.67  cm/s 

1.39 

1.46 

1 . 76  cm/s 

Oscar 

*  1.52  an/s 

1.42 

1.42 

1 . 62  rm/s 

BEC 

*  75  cm 

69 

175 

133  cm 

l.ageos 

*  8  cm 

8 

8 

16  cm 

Lageos 

*  7  cm 

7 

7 

17  cm 

Geos-2 

*  134  cm 

85 

97 

347  cm 

*  see  text 


The  standard  field  used  in  Table  26  is  a  specific  field  that  in  some  cases 
is  a  tailored  geopotential  model.  Specifically  the  standard  models  are  PGS1331 
for  Starlette;  PGSS4  for  Seasat;  GEM10B  for  Oscar,  BEC,  and  Geos-2;  GKML2  for 
l.ageos.  The  86C/D  and  GPM2  models  were  used  in  the  fits  complete  to  degree 
and  order  36. 


9.  Summary  and  Conclusions 

It  has  been  five  years  since  our  previous  high  degree  expansion  and 
combination  solution  had  been  carried  out.  Since  then  we  have  seen  a  number 
of  developments  that  warranted  the  new  solutions  described  in  this  report. 
These  developments  are  in  the  data  area;  the  theoretical  area;  and  the 
computer  area. 

In  the  data  area  wo  have  recently  completed  a  new  set  of  l'xl'  mean 
terrestrial  gravity  anomalies.  This  collection  of  anomalies  has  increased  from 
our  prior  set  and  has  improved  estimates  in  a  number  of  areas.  In  addition  a 
set  of  ocean  wide  l'xl*  anomalies  has  been  derived  from  the  Geos-3/Seasat 
altimeter  data.  The  merged  terrestrial/altimeter  data  set  provides  a  more 
complete  and  more  reliable  data  had  than  used  before.  And  finally,  we  had 
the  GEML.2  potential  coefficient  set  and  its  revised  accuracy  estimates  to  use 
for  our  a  priori  solution. 

In  the  theoretical  area,  several  improvements  have  been  made  in  the  new 
solution.  First,  the  boundary  condition  relating  gravity  anomalies  and  the 
disturbing  potential  was  more  precisely  formulated  to  avoid  a  spherical 
approximation.  Second,  correction  terms  were  formulated  that  enabled  the 
integration  over  the  surface  of  the  ellipsoid  as  opposed  to  an  integration  over 
a  spherical  surface. 


3.  What  is  the  smoothing  effect  on  the  potential  spectrum  from  the  use  of 
the  optimal  estimation  (or  least  squares  collocation  procedure)? 

4.  What  is  the  best  way  to  carry  out  a  high  degree  expansion  using 

30  x30  mean  values  and  what  differences  will  be  found  with  the  1*  solution 

carried  out  for  this  study0 

5.  Sould  Doppler  derived  geoid  undulations  be  incorporated  into  the 
geopotential  solution  and  if  so,  what  parameters  should  be  modeled? 

6.  What  is  the  best  way  to  balance  the  needs  for  a  gravity  near  the 

earth’s  surface  and  in  space? 

7.  Should  geophysical  anomalies  be  used  at  all  in  the  solutions  or  should 

their  use  be  restricted  to  locations  where  some  independent  verification  is 

possible0 

8.  What  is  the  reason  for  the  residual  undulation  correlation  with 
elevation  found  in  most  geopotential  solutions? 

9.  What  is  the  effect  of  using  the  full  error  variance-covariance  matrix  of 
the  a  priori  potential  coefficients. 

10.  What  are  the  divergence  problems  for  the  potential  series,  and  series 
representation  of  other  gravimetric  quantities,  implied  by  these  high  degree 
spherical  harmonic:  expansions0 

Finally  we  should  note  that  there  is  a  strong  need  for  an  improved 
satellite  potential  coefficient  set.  Although  GEML2  is  strong  at  degrees  below 
6,  substantial  improvement  is  needed  in  the  higher  degrees.  As  such  fields 
become  available  it  will  become  reasonable  to  repeat  the  combination  process  to 
obtain  an  accurate  representation  of  the  earth’s  gravitational  potential. 
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